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Abstract. We study some new universal aspects of diffusion in chaotic sys¬ 
tems, especially such having very large Lyapunov coefficients on the chaotic 
(indecomposable, topologically transitive) component. We do this by dis¬ 
cretizing the chaotic component on the Surface-of-Section in a (large) number 
N of simplectically equally big cells (in the sense of equal relative invariant 
ergodic measure, normalized so that the total measure of the chaotic compo¬ 
nent is unity). By iterating the transition of the chaotic orbit through SOS, 
where j counts the number of iteration (discrete time), and assuming com¬ 
plete lack of correlations even between consecutive crossings (which can be 
justified due to the very large Lyapunov exponents), we show the universal 
approach of the relative measure of the occupied cells, denoted by p{j), to 
the asymptotic value of unity, in the following way: p{j) = 1 — (1 — A)i^ so 
that in the limit of big N, N ^ oo, we have, for j/N fixed, the exponential 
law p{j) 1 — exp{—j/N). This analytic result is verihed numerically in a 
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variety of specific systems: For a plane billiard (Robnik 1983, A = 0.375), for 
a 3-D billiard (Prosen 1997, a = —1/5, b = —12/5), for ergodic logistic map 
(tent map), for standard map {k = 400) and for hydrogen atom in strong 
magnetic field (e = —0.05) the agreement is almost perfect (except, in the 
latter two systems, for some long-time deviations on very small scale), but 
for Henon-Heiles system [E = 1/6) and for the standard map {k = 3) the 
deviations are noticed although they are not very big (only about 1%). We 
have tested the random number generators (Press et al 1986), and confirmed 
that some are almost perfect (ranO and ran3), whilst two of them (rani and 
ran2) exhibit big deviations (in units of the theoretically expected standard 
deviation, as much as 20). Thus physical deterministic chaotic systems can 
be better simulators of random numbers than some well known mathematical 
algorithms. We give an outline of an improved analytical theoretical model 
(the so-called two and many component model), where deviations from the 
exponential law can be captured in a statistical way. 


PACS numbers: 05.45.-|-b, 05.40.-|-j, 05.60.-l-w, 03.20.-|-i 
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One of the major open problems in the mathematics of nonlinear Hamil¬ 
tonian (conservative) chaotic systems of the KAM type is the proof of the 
so-called coexistence problem (Strelcyn 1991), i.e. the proof that the chaotic 
components have positive measure. (The KAM Theorem guarantees that the 
set of invariant tori has positive measure, whose complement is small with 
the perturbation parameter (Kolmogoroff 1954, Arnold 1963, Moser 1962, 
Benettin et al 1984, Gutzwiller 1990).) The chaotic component could be de- 
hned e.g. by the positivity of the (largest) Lyapunov exponent, which is a 
sufficient but not a necessary criterion]^. We shall define a chaotic component 
as the closure of a dense chaotic orbit, which is thus assumed to be an inde¬ 
composable invariant component (topologically transitive, i.e. containing a 
chaotic dense orbit). 

There are really no serious doubts about the positivity of the measure of 
the chaotic components, and so in physics we rely on heuristic arguments to 
actually assume positivity. Then the question is how to calculate the sym- 
plectic (invariant and ergodic) measure of the chaotic component. 

We have approached this problem in a recent extensive work (Dobnikar 1996, 
Robnik and Dobnikar 1997) on the dynamics in a plane billiard system, de¬ 
fined as the quadratic conformal map of the unit disc (in the complex z-plane) 
onto the physical (complex) tc-plane, w{z) = z + Xz"^. This system has been 
introduced by Robnik (1983), and recently extensively studied for many dif¬ 
ferent values of A by many authors, in a variety of contexts and even in 
experimental setups like quantum dots (Bruns and Stone 1994, Stone and 
Bruns 1993ab), optical model (Nockel and Stone 1997, Nockel et al 1996), 
microwave cavities (Rehfeld et al 1996, Richter 1996, Stockmann et al 1997). 
Further dynamical details were corrected in (Hayli et al 1987). 

At A = 0 we have the integrable case of the circle billiard, for 0 < A < 1/4 
the billiard is convex, and since the boundary is analytic, the KAM theory 
applies (Lazutkin 1981, 1991), at A = 1/4 we get the first point of zero 
curvature at the boundary point w = w{z = —1) = —3/4, allowing for the 
breaking of Lazutkin caustics (invariant tori associated with the boundary 

^For example, in nonrational plane polygonal billiards all Lyapunov exponents are 
strictly zero (Sinai 1976), and yet they can be ergodic. 
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glancing orbits), and for 1/4 < A < 1/2 the billiard is non-convex, largely 
and strongly chaotic (very tiny islands of stability), probably becomes rig¬ 
orously ergodic at some A > 0.2775 (Li and Robnik 1994), and is definitely 
proven to be rigorously ergodic for A = 1/2 (the so-called cardioid billiard, 
having the cusp singularity aX w = w{z = —1) = —1/2) (Markarian 1993). 
The cardioid billiard has been studied also by Backer, Steiner and Stifter 
(1995) and by Backer and Dullin (1997). 

Our main problem was to calculate numerically, accurately and reliably, the 
measure of chaotic components. Working in the KAM regime (0 < A < 1/4) 
we were observing the typical KAM hierarchy of smaller and smaller islands 
of stability surrounded by chaotic components, the details of which will be 
reported in a separate paper (Robnik and Dobnikar 1997), but are reported 
already in (Robnik 1983). 

The main objective was to calculate the fractional measure (the relative 
area on the SOS, the latter being dehned by the PoincarABirkhoff coordi¬ 
nates) of the largest chaotic component at given A, which we traditionally 
denote by p 2 - 

This parameter is important also in treating the related quantum mechan¬ 
ical problem (solutions of the Helmholtz equation with the Dirichlet bound¬ 
ary conditions on the billiard boundary) (Berry and Robnik 1984, Prosen 
and Robnik 1993,1994a-b, Li and Robnik 1995). 

We have discovered, to our big surprise, that this numerical calculation is 
extremely difficult, and as one consequence, some of the previous results had 
to be revised. By dividing the SOS into a large number of rectangular grid 
cells of equal relative (normalized) measure we have calculated the relative 
measure of the chaotic component p 2 by three different methods: (Ml) 
calculating the Lyapunov exponent for a trajectory starting in a given cell and 
summing up the area of cells having the positive Lyapunov exponent; (M2) 
calculating two nearby trajectories, separated by an inhnitesimal distance 
(e.g. single precision e-8 while all calculations were in double precision e- 
16) and summing the cells exhibiting macroscopic divergence in a reasonable 
time; and (M3) starting the chaotic trajectory and counting (summing up) 
the area of the cells visited (’’black cells”). 

It turned out that the third method (M3) is the best, fastest, and most 
reliable. To improve the result one has to enlarge the number of cells N 
on the chaotic component. But then, the time j (of iterating the map on 
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the SOS, the discrete time) has to be taken much larger (by at least several 
orders of magnitude) than N, otherwise the statistics of visiting cells (black 
cells) would be insignihcant. Even after calculating as many as 10® itera¬ 
tions, with N = 1000 — 2000, the result was no better than within 1%. One 
reason is that the boundary of the chaotic component turned out to have 
relatively large fractal dimension around 1.56. The difficulties of estimating 
the asymptotic value of the relative area p 2 of the chaotic component led 
us to the careful investigation of the evolution of relative area of occupied 
(black) cells P 2 (j) as time j proceeds, and we wanted to understand that 
theoretically, in order to be able to make better estimates of the limiting 
asymptotic value p 2 = P 2 U = cxd). From here onwards we shall use the no¬ 
tation p{j) = p 2 {j)/p 2 {j = 00 ). 

In general this time evolution p{j) with j is quite complex and specihc, 
nonuniversal, depending on many features appearing in the phase space 
(SOS), e.g. existence of sticky objects like cantori can affect a temporary 
but quite persistent trapping of the orbit near such an object, which is then 
manifested in a transient plateau of the curve p{j), which sometimes might 
be mistakenly interpreted as final and definite convergence of the cumula¬ 
tive area/volume P 2 {j)- However, if the system is really strongly chaotic, 
having large maximal positive Lyapunov exponent, then due to the bound 
motion and conservation of the phase space volume, we hud a very strong 
stretching and folding in the phase space. In such a limiting case therefore 
one small phase space cell (SOS cell) becomes uniformly distributed, in the 
coarse grained sense, all over the allowed chaotic component. Thus, in such 
extreme case, the probability of entering a given cell belonging to the same 
chaotic component in SOS is just equal to the relative measure of the cell, i.e. 
there are no correlations, not even between two consecutive SOS iterations: 
complete randomness of deterministic motion. 

While the behaviour in such ideal extreme case is quite obvious, it is far 
from obvious that the conditions of the complete randomness are actually 
satisfied in specihc chaotic deterministic dynamical systems. Therefore, to 
make things precise we have developed the following theoretical model. 

Suppose we have N cells, where their order and geometry of arrangement 
is completely irrelevant and perhaps not even known. We are hlling the cells 
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with balls, one at the time. At each step j (discrete time) we have equal 
probability to choose any cell, equal to a = 1/N. (Thus there are abso¬ 
lutely no correlations between any moves, including the consecutive ones, 
and therefore e.g. the repetition of falling into a given already hlled cell is 
allowed.) We dehne by Pj{k) the probability that at j-th step k cells are 
occupied, keeping in mind that a = 1/N is the model parameter implicit 
in the mathematical formulae. We shall refer to this model as the random 
model (of strongly chaotic deterministic diffusion). The probabilities must 
be normalized, therefore 


k=N 

Y,P,{k) = l, Vj = l,2,.... (1) 

fc=l 

We shall calculate Pj{k), and their moments, in particular the hrst moment, 
namely the average (normalized) measure of the occupied cells, 

N 

pU) = H kaPj{k) = (ka), (2) 

k=l 

where by (...) we denote the averaging operation. 

Before explicitly calculating Pj{k) we observe the physically (probabilisti¬ 
cally) quite obvious recursion relation, namely 


P,+i(«: + l) = P,(A^+l)f^ + P,W(l-k), V0<A^<j, (3) 

where we also dehne the boundary conditions 

P,(0) = 0, Pi(l) = l, P,{k>j) = 0, Pj{k>N) = 0. (4) 

The interpretation of equation (^ is: the probability to have {k + 1) cells 
occupied at time (j -|- 1) is equal to the sum of the following probabilities: 
either at time j the {k + 1) cells were already occupied, and we add the 
next ball into the black cells with probability {k + 1)/N, or at time j only 
k cells are occupied (black), and we add the next ball (the (j -|- l)-st one) 
into the empty (not-yet-occupied) cells with probability (1 — k/N). With 
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the boundary conditions the recursion relation (^) solves the problem, 
in principle. We show the explicit solution below, using different approach. 
However, for a practical (numerical) evaluation of Pj{k) (on the computer) it 
is much better to use the recursion formula (j^) than the explicit result which 
we shall derive below. 

For the beginning please observe that the summation of the recursion equa¬ 
tion d^) on each side from k = Q io k = N — 1 confirms the preservation of 
normalization (|lD, for all j. 

Next, we can find the solution for p{j) at once by the following trick: multi¬ 
ply the recursion relation (|^) on the left and on the right by (A; -|- l)/iV and 
sum it from A: = 0tofc = iV — 1. By denoting S{j) the second moment, 

k=N 

S(f) = {(ak-f) = x: P,(k)(ak)\ (6) 

k=l 

we obtain in a straightforward manner. 


p(j + l) = SU)-{SU)-P,{N)) + {l-a){pU)-P,{N))+a{l-P,{N)), (6) 

and therefore after cancellation of S{j)'s we get the simple recursion equation 
for p{j), namely 


p{j + l) = a + {l-a)p{j), (7) 

with the explicit solution, quite easy to find, 

P(i) = 1-(1-a)^ = 1-(1-(8) 

which in the limit of large N, for j/N fixed, becomes simple exponential law 

p(j) ~ 1-exp(-^). (9) 

We see that in the limit of sufficiently large N, for j/N fixed, we have the 
universal scaling of the relative measure of chaotic region p, normalized to 
unity, such that p{j) —>■ 1, when j —>■ oo, as a. function of the scaled discrete 
time, namely j/N. We shall show and see below, that this law is obeyed by 
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a surprising variety of deterministic dynamical systems. 


Our random model is a probabilistic (statistical) model and therefore we 
can calculate all moments of Pj{k), systematically, using the same trick as 
above: multiply the recursion relation by {k + on both sides, sum 

it up from fc = 0tofc = A^ — Ion both sides, and uncover the recursion 
relation for the second moment S{j), namely 




( 10 ) 


and by using the exact result for p(j) from equation we have 

5(j + l) = 2a-a(2-a)(l-a)^' + (l-2a)^(j), a = l/N. (11) 

This equation can be solved either by standard technique or by using the 
dehnition of S{j), equation (^, to yield the explicit result 


^(j) = l-(2-a)(l-a)^ + (l-a)(l-2a)^ a = l/N, (12) 


so that the predicted dispersion cr^(j) according to our model is exactly 


^■^(j) = >S'(j)-p^(j) = a(l-a)^ + (l-a)(l-2a)'^-(l-a)^^ a = 1/iV. (13) 


In the asymptotic limit of sufficiently large number of cells N = 1/a —^ oo, 
but keeping j/N hxed, we hnd the simple exponential laws: 


p{j) 1 - exp(-j/A^), S{j) ^ [1 - exp{-j/N)]‘^ ^ (14) 


and therefore 


K N '[exp(-VA') - exp(-2j//V)| ^ 0. 


(15) 


Now we show the explicit and exact result for Pj{k), just for the sake of 
completeness. In fact the quantity we seek is the subject of the classical 


problem from combinatorial analysis, treated e.g. in (Vinogradov 1979, Vol.2, 
p.973, Riordan 1978, p.48. Table 2, Graham, Knuth and Patashnik 1994): 
The question is how many possibilities are there to distribute j different 
things into N different cells under the condition that N — k cells are empty 
(i.e. precisely k cells are occupied): The answer is well known, namely in the 
literature denoted by Cp^j{N — k), 

C^,(N-k)=(^j^^_^yiS(3,k) = j^;j^SU,k), ( 16 ) 

where S{j,k) are the so-called Stirling numbers of the second kind (Vino¬ 
gradov 1979, Riordan 1978, Graham, Knuth and Patashnik 1994) 

(IT) 

which are known to satisfy the triangular recursion relation S{j, k) = kS{j — 
1, A;)-|-iS(j — l,/c — l), where iS(0, 0) = 1, 5(j, 0) = 0, iS(0,A;) = 0 {j,k)> 

0. Of course, having the N cells and j things, like in our random model, the 
total number of possibilities to distribute j things (balls) into N cells is just 
, and therefore we have the hnal and complete explicit solution to the 
random model, namely 


^ c^,{N-k) ^ msu, k) 

^ Ni {N-k)\N^' 


(18) 


Now we proceed by analyzing specihc dynamical systems from the point of 
view of the statistical theory presented in our random model, to see to what 
extent we hnd agreement in real systems. The really big surprise is that 
the behaviour was found in excellent or even perfect agreement with theory 
in a large variety of deterministic dynamical systems, sufficiently far from 
a pronounced KAM-regime, by which we mean either close to ergodic (big 
P 2 (j = cxd) 1), or strongly chaotic (big Lyapunov exponent but not neces¬ 
sarily very large p 2 (j = oo)). 

We have investigated the plots p{j) versus j/N for the following systems: 
(a) the 2-D billiard (Robnik 1983) at A = 0.375, (b) the 3-D billiard (Prosen 
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1997, def. geom.: a = —1/5,6 = —12/5), (c) the hydrogen atom in strong 
magnetic held (DKP = diamagnetic Kepler problem) with the scaled energy 
e = —0.05 (see e.g. Hasegawa, Robnik and Wunner 1989), (d) the Henon- 
Heiles system at the escape (dissociation) energy E = 1/Q (see Henon and 
Heiles 1964), (e) the standard map (Chirikov 1979) with k = 3, (f) the 
standard map (Chirikov 1979) with k = 400, and (g) for the logistic map at 
A = 4 (the ergodic tent map). 

In case of smooth systems (c) and (d) we have used special symplectic 
integration routines, devised by Yoshida (1990). This enabled a fast and 
extremely accurate calculations, allowing us to compute about the order of 
magnitude 10^ iterations on the SOS. 

The agreement for all systems was perfect (deviations much smaller than 
1%), except in (d) and (e), where the deviations huctuated around up to 
1%. Therefore we do not show these plots, since all curves practically overlap 
with the theoretical curve (H) and (P) within the graphical resolution. The 
small deviations stem from the fact that the Lyapunov coefficients are still 
not big enough, and also that there might be signihcant episodes of tran¬ 
sient behaviour in the relationship of p{j) with respect to the discrete time 
j/N. Such transient episodes are typically caused by the sticky objects in the 
phase space, e.g. by cantori, where the classical orbit spends long time be¬ 
fore resuming the chaotic random hlling of the remaining empty cells. They 
are very well manifested in systems with more pronounced KAM structure, 
like e.g. the 2-D billiard (A = 0.15) or 3-D billiard (for sufficiently small a 
and b, such as a = —0.1, 6 = 0). In order to describe such systematic ef¬ 
fects in a statistical way we have developed a multicomponent random model, 
where orbital transitions inside each component are random as in our random 
model, however, they might jump (rarely) from one into another component. 
Some of the dynamical features are well described by such a model, whose 
detailed description will be published in a separate paper (Robnik, Prosen 
and Dobnikar 1997). 

In cases (a) - (c) the agreement is so good, quite surprisingly, that it is nec¬ 
essary to magnify the scale so that the details of deviations become clearly 
visible. This is done in hgures l(a-c) for the systems, respectively. We 
plot the fluctuation (difference) PnumericaiU) - PtheoryU) (the noisy curves) 
to be compared with the theoretically expected standard deviation a{j) = 
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^\/S{j) — p^{j) (smooth curves), again as a function of the scaled discrete 
time j/N. We do this in each of the plots for one initial condition and for 
the average over 50 randomly (uniformly over the chaotic component) chosen 
initial conditions which suppress the dispersion by a factor of 50, of course, 
and the standard deviation by 

The conclusion in inspecting these plots is that the 2-D billiard perfectly 
obeys the law of our random model, whilst for the hydrogen atom in a strong 
magnetic held and the 3-D billiard we uncover systematic deviations from 
the theory for sufficiently large scaled discrete times j/N: the deviations are 
orders of magnitude bigger than the prediction of our random model for the 
DKP, but somewhat smaller in the 3-D billiard. It should be noted that the 
deviations are almost strictly negative, rehecting the fact that the physical 
orbits like to stick to already occupied cells (existence of sticky objects in 
the phase space). 

It is interesting to look at the results for random number generators. Cells 
are ’’randomly” chosen, using different random number generator algorithms. 
One such algorithm was devised by Finocchiaro et al (1993) in the context 
of nuclear physics. The agreement was perfect and so we do not show the 
fluctuations plots. 

Also, we have checked and tested some well known algorithms for random 
number generators, such as ranO, rani, ran2 and ran3 devised and described 
in (Press et al 1986, ch. 7, and the references therein). The results are 
shown in hgure 2. As we see ranO and ran3 are in excellent agreement with 
our random model, whilst for rani and ran2 random generators the deviations 
become clearly very large. Thus, for example, our dynamical deterministic 
2-D billiard systems, or the logistic map, or the standard map at k = 400, 
are in fact better random number generators than some of the most familiar 
random number generators used in the computers. 

Finally, we have looked at the irrational triangle (the angles are a = {\/5 — 
l)7r/4, (3 = (v^ — l)n/2) as one ergodic system with strictly zero Lyapunov 
exponents (Sinai 1976). Even there the agreement with our random model 
is surpringly good on the largest scale (hgure 3), while the huctuation di- 
agramme exhibits large deviations, again negative, showing that the real 
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billiard orbits like to stay in the already occupied regions (sticky objects in 
phase space). 

In conclusion, we have developed the random model of stochastic diffusion of 
dynamical systems with invariant measure on their Surfaces of Section, which 
is supposed and conhrmed to apply very well in strongly chaotic systems (for 
which the Lyapunov coefficients on the chaotic component are sufficiently 
large). We have discovered and explained the universal scaling behaviour of 
the normalized chaotic measure p as a function of the scaled discrete time 
j/N, where N is the number of cells: namely, in the limit of sufficiently large 
N, ioT j/N fixed, we have a simple exponential law p{j) = 1—exp(—j/iV). We 
also predict the higher moments, especially the dispersion given in equation 
(|T^. The model is solved completely in the sense that we have calculated the 
probabilities Pj{k) of having exactly k non-empty cells at time j. Therefore, 
one can calculate all the moments. The deviations from the predictions of 
the random model are qualitatively understood, but will be treated in detail 
together with a new, more general theory (the multicomponent model) in a 
separate work (Robnik, Prosen and Dobnikar 1997). The work is in a sense 
extension of the theory of transport in Hamiltonian systems by MacKay, 
Meiss and Percival (1984). 
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Figure captions 


Figure l(a-c): We show the plots of PnumericaiU) “ PtheoryU) versus the 
scaled discrete time j/N. The outer noisy curve is the numerical result for a 
chaotic orbit with a certain representative initial condition, whilst the inner 
one is the average over hfty evenly distributed initial conditions. We also 
show the ±cr(j) standard deviation, as predicted theoretically in equation 
(13). It is clearly seen that for the 2-D billiard in (a) the agreement is almost 
perfect, in the sense that the fluctuations are within the predicted range, 
whilst for the 3-D billiard in (b) we see systematic deviations, obviously 
caused by some perhaps unexpected long time correlations. Snch correlations 
are even stronger in case of the hydrogen atom in a strong magnetic held, 
shown in (c). In all cases the deviations are predominantly negative: the 
orbit tries to stick to some of the already occupied cells. 


Figure 2; We show the results for the random nnmber generators, compared 
by the theoretical ±cr(j) curves according to equation (13) (Press et al 1986): 
for ranO and ran3 the agreement is excellent, while for rani and ran2 the 
deviations are so big, that they are actually not random: they seem to repel 
from the occupied cells. 


Figure 3: We show the global plot, Pnumericai{j) (full line) and ptheoryU) 
(dashed), for the irrational triangle, where the agreement is snrprinsingly 
good, in spite of strictly vanishing Lyapnnov exponents. In the huctuations 
diagramme we plot Pnumericai{j) - PtheoryU) versus j/N, together with the 
±(j, where we see the same trend to negative deviations dne to the sticky 
objects in the phase space. 
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